function f = fSpring(p1,p2,k,l0) 
% function f = fSpring(p1,p2,k,l0)
% Spring Force for Linear spring
% k is stiffness of each kind of spring
% l0 is intrinsic length of all springs (vectorised)
% p1 and p2 are vectorised set of points specifying endpoints of springs
% p1,p2 and l0 are all vectors of same lengths

f = zeros(size(p1)); %  [Magnitude Direction]
dx=p2-p1;
dist=sqrt(sumsq(dx'));
fTheta=angle(complex(dx(:,1),dx(:,2)))';
fMag = - k* (dist - l0) ;
f(:,1) = fMag.*cos (fTheta);
f(:,2) = fMag.*sin (fTheta);
